Emergence of coherence in the Mott-superfiuid quench of the Bose-Hubbard model 
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We study the quench from the Mott to the superfiuid phase in the Bose-Hubbard model and 
investigate the spatial-temporal growth of phase coherence, i.e., phase locking between initially 
uncorrelated sites. To this end, we establish a hierarchy of correlations via a controlled expansion 
into inverse powers of the coordination number 1/Z. It turns out that the off-diagonal long-range 
order spreads with a constant propagation speed, forming local condensate patches, whereas the 
phase correlator follows a diffusion-like growth rate. 
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A sweep through a symmetry-breaking quantum phase 
transition is one of the simplest ways to create and am- 
plify quantum correlations, i.e., entanglement. As the 
initial quantum state is symmetric, all directions of sym- 
metry breaking are equally likely and seeded by quan- 
tum fluctuations. Furthermore, the diverging response 
time at the critical point indicates that the many-particle 
quantum system is driven far away from equilibrium dur- 
ing the sweep. While nearby points will most likely break 
the initial symmetry in the same direction, two very dis- 
tant points may spontaneously select different directions 
of symmetry breaking As a result, the spatial or- 
der parameter distribution after the quench will be in- 
homogeneous - and its spatial correlations are directly 
determined by quantum correlations |2|]. 

In contrast to the static properties of classical and 
quantum phase transitions, much less is known about 
their non-equilibrium dynamics. Motivated by recent ex- 
periments [3, 4] , as well as from a more fundamental point 
of view, this field of research is attracting increasing in- 
terest [5(. Open questions in this context include: How is 
the new order parameter established and how fast does 
it spread? Are there universal scaling laws [6] similar 
to those in static phase transitions close to the critical 
point? In this Letter, we address these questions for the 
Bose-Hubbard model which is considered Q one of the 
prototypical examples for quantum phase transitions and 
also relevant for experiments in optical lattices The 
Bose-Hubbard Hamiltonian is given by (fi, = 1) 



J 



11 = -^zZ T ^l^ 



U 
~2 



(1) 



where aj^ and d„ are the creation and annihilation opera- 
tors at the lattice sites p and v, respectively. The lattice 
structure is encoded in the tunneling matrix T^„ € {0, 1} 
and J denotes the hopping rate. The coordination num- 
ber Z = T^u € N counts the number of tunnelling 
neighbors at any given site p. Finally, U is the on-site 
interaction with — a^a^. For simplicity, we assume 
an average filling of one boson per site (n^) = 1. 

At a critical ratio of J/U (see below), this model ([T]) 
features a quantum phase transition separating the sym- 



metric Mott insulator phase with a gap from the super- 
fluid phase with long-range order (i.e., broken symme- 
try) and Goldstonc modes [8| . We start deep in the Mott 
phase J = with |\&Mott) = Yifi Op |0) an d quench to the 
superfiuid regime by means of an instantaneous switch- 
ing to a finite value of J. The theoretical description of 
the ensuing non-equilibrium many-particle dynamics as- 
sociated with this process is a non-trivial task for the 
non-integrable Bose-Hubbard Hamiltonian (JlJ. There 
are two major options: numerical computations which 
are limited to systems of finite size or analytical calcu- 
lations which require suitable approximations. In order 
to control the accuracy and consistency of these approxi- 
mations, they should be based on the expansion in terms 
of a small or large parameter. For the Bose-Hubbard 
model, this could be a large filling (n M ) 3> I [SJ or a large 
number of interaction partners for the extended version 
In this Letter, we assume a large coordination num- 
ber Z> 1 and employ a systematic expansion in 1/Z. A 
large Z can occur in a large number of spatial dimensions 
or a large number of tunnelling partners. In the follow- 
ing, we focus on the second case and assume two spatial 
dimensions for simplicity. 

Let us consider the reduced density matrices for one 
lattice site p^ = Tr^{p} and for two sites p Mt , = Tr^{p} 
etc. Furthermore, we separate the correlated parts via 

Pn» = P\v + PvPv as wel1 as Pl^\ = P% X + PUvfa + 

Pli\Pv + PuxPn + P^Pvf>\ etc. Our derivation is based on 
the following scaling hierarchy of correlations 



Ps 



q( Z 1-\S\ 



(2) 



where |<S| is the number of lattice sites in the set <S, i.e., 
p, = 0(Z°), pl v = 0{Z-% f>l vX = O(Z^) etc. This 
hierarchy ([2]) is similar to the quantum de Finetti the- 
orem [llj], the generalized cumulant expansion [12j, and 
the BBGKY hierarchy but we are considering lattice 
sites instead of particles. In order to derive the hierarchy 
([2|), we introduce the generating functional 



T{a) = F({a^}) = In Tr 




(3) 
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where d M is an arbitrary operator acting on the lattice 
site p,. This functional generates all correlated density 
matrices via d M -derivatives ft = dT/da^l^o, as well as 
Puv = d 2 F / da^da^a^Q etc., where we have used the 
notation (n^\ dF/da^ \m^) — dT/d (n M | d M \mftj. In- 



troducing the Liouville super-operators £ M and L^ LV via 

id t p = [H,p] = £ M + £^p/ z , the temporal 

evolution of J 7 is given by 



id t F{&) = Tr M |d M ^J^} + 7? E Tl > { 



(a M + a v + d M d y ) 



/ d 2 T dT 8T \ 
\ daada u dan <9dv / 



(4) 



By taking successive derivatives, we establish the following set of equations for the correlated density matrices 

1 l r ■pup=5\{ M } 



lies 

VvjP=S\{ti.,v} 



}U7> 



VCS\{n} 
QUQ=P 

Qcv 



(5) 



M,i>S<S VCS\{fi,iy} 

where C s — + C vtl . A careful inspection of this set of equations shows that the hierarchy in © is preserved 



in time. If all the correlated density matrixes ftp on the r.h.s. obey the hierarchy ftp = 0(Z 1 ~l' p l), then the time 
derivative on the l.h.s. does also satisfy ©. Therefore, starting deep in the Mott phase I^Mott) = 11^ ®ti |0) where © 
is trivially satisfied since all correlations vanish, we find that © remains valid for a finite time; more precisely, for a 
time scale of 0(lnZ) limited by the instability of the growing modes, see below. Let us consider some examples for 
the above evolution equation: for one lattice site S — {/i}, we get 



id t ft = C^ft + ^ E Tl ' K { £%(pU + PvP*)} 



(6) 



Since ft K is suppressed by 0(1/Z), we may neglect this term to the lowest order in 1/Z and thereby obtain the 
Gutzwiller approach [l4|, EH- Starting in the Mott phase (a re ) = 0, we have Tr K {£^ K p K } = and thus we get 
ft = |1) (1| + 0(1/Z) = ft + 0{1/Z). The two-point correlation can be studied with S = {p, v} 



idtft^ = C^plv + ^£^(ft^ + ftft) + E TrK {£%(P»™ + P%P* +P C ukPiS\ ~ ^ Tr ^ \^l^P% + PuM] 



z 

+(p O v) . 



(7) 



Again, to the lowest order in 1/Z, we may neglect the three-point correlation ft UK = 0(1/ Z 2 ) and replace ft by p° 
arriving at a closed set of equations for ft . Furthermore, using Tr„{ C^ft} — 0, this simplifies to 



id t pl v = ( 4 + £,) P% + ~ CUIpI + I E Tr « { KnPlnt + ZlnPlJl) + 0{1/Z 2 ) . 



(8) 



r 



|0) M (1| as local 



Introducing p^ — \1) ^(2\ and 

particle and hole operators, we find that their correla- 
tion functions = Tt{ph\K}, = Tr{p 

III = ^{pP^hu}, and f™ = Tr{pp^p v }, obey a closed 
linear system of equations. Correlation functions 
containing higher occupation numbers n, m > 3 decou- 
ple and obey a homogeneous set of equations without 
the source terms stemming from C^ft^ft. Thus, they 



are trivially zero assuming the Mott state initially - up 
to the accuracy 0(1/Z 2 ) under consideration. 

Assuming discrete translational invariance, we may 
simplify the equations via a Fourier transformation 

(id t -U + 3JT k )fl 2 = -yftmup + ff + i), 

idtfl 1 = *d t f 22 = V2JT k (fl 2 - fl 1 ) , (9) 



With /: 



21 



12\* 



We find that f^ 1 = /| 2 indicating an 



3 



effective particle-hole symmetry. The evolution is fully 
determined by J, U and the Fourier transform of the 
tunnelling matrix 



+ 0(k i ). (10) 



Assuming discrete rotational symmetry of the lattice, we 
obtain exact isotropy at small k with a unique effective 
mass m*; otherwise one would have m* =^ m* etc. For 
large Z, the sum over lattice sites /x in (|10|) involves many 
terms and thus only small wavenumbers yield significant 
contributions, corresponding to the effective mass being 
small m* = 0(1/ Z). 

The linear system in ([9]) can be rewritten in matrix 
form idtfk = [£fc,ffc] + Sfc with f fe = (/£ m ), the source 
term Sk and the Liouvillian matrix , whose eigenvalues 
uj^ h are associated with quasi-particle/quasi-hole exci- 
tations [3]. As one would expect from the particle- hole 
symmetry, two of the four eigenfrequencies of the full 
linear system in ([9]) vanish and the other two read 



U 2 - QJUT k 



J 2 n 



(11) 



After the quench from J- m <C U to J out = J, the set 
of equations in © is solved for the initial conditions 
(hl~hv)o = S^u and (frjjp„) = {P\K)o = (p^)o = 0. 



Inserting a p 



V2 



where the remaining 



terms containing higher occupation numbers do not 
contribute, we obtain 



a\(t)a v (t)) 



ST ik-(Tu-Tu) 4JUT k 1 - COS^fet) 

V N "2 



where AT denotes the total number of lattice sites. 
This prediction could be experimentally verified. After 
preparing bosonic atoms in an optical lattice deep in the 
Mott state J = and quenching it to a finite J ou t = J, 
time-of-flight images for varying time intervals t after the 
quench yield the Fourier transform of (cl 1 (t)a v (t)) and 
allow the determination of uik In the Mott regime 
>^out < Jen this sweeping method allows not only the 
measurement of the gap value Wfe = o which is a signature 
of the insulator phase, but also the k dependence of the 
energy uik required for particle-hole pair creation. 

The expression under the square root in (fTTj) has two 
zeros J± = {7(3 ± \/8) for k — 0. The first one marks 
the critical point J cr = J_, see, e.g., [14]. For J < J cr , 
i.e., in the Mott phase, all modes are stable Uk 6 K. 
For J CT < J < U (which is in the superfluid regime), 
all modes with wavenumbers below the value fc cr given 
by JX/c cr = J cr become unstable and the mode k = 
yields the fastest growth. For J > J7, on the other hand, 
a finite wave-number fc* > yields the fastest growth 
and dominates the evolution of (fit (^fi^i)), see Eq. (JTSJ) 
below and 16]. Finally, for J > J + , long- wavelength 



modes become stable again Wfc = o € K and only the modes 
within a finite fc-interval grow. 

In order to search for universal behavior close to the 
critical point, we study a quench not too far into the 
superfluid regime, i.e., J = J CI (1 + e) with < e -C 1. 
In this case, the dispersion curve lu 2 , dips below zero for 
small k only, and thus the growth of correlations can be 
determined using the long- wavelength approximation 



2 k' 



(13) 



where c 2 = 3J(U — J)/m* = O(Z) is a velocity scale. 
Note that the growth rate j 2 ~ J — J cr strongly depends 
on the distance to the critical point, whereas c 2 is nearly 
constant. For large lattices N ^> 1, the sum over fe is 
well approximated by an integral. For large jt ^> 1, 
this integral becomes dominated by the fastest growing 
modes and, thus, can be estimated using the saddle-point 
approximation 

(fit (t)6 M (*)) « Af(t) exp | 7v /t 2 - - tv) 2 /c 2 j ,(14) 

where the time dependence of the normalization factor 
Af(i) — 0(1/ Z) is weak (power-law) compared with the 
exponential growth in (|14[) . Focusing on the dominant 
exponential part in (|14p . we find a constant propagation 
speed c of the correlations similar to the Licb-Robinson 
bound 17[ . Furthermore, there is a universal scaling be- 
havior. Moving towards or away from the critical point 
through the rescaling 7 — > 7' = A7, an analogous rescal- 
ing of time t — > t' = t/X and distance r —> r' = r/X 
leaves the dominant behavior of (|14p invariant. 

As another application of (fl"4"|) . let us determine the 
condensate fraction within a compact lattice region <S 
of size \S\ ^> 1. By analogy with the continuum case 
fjTJ . the condensate fraction is defined via the largest 
eigenvalue of the two-point correlation function (a*Ji v ). 
Within S, the largest eigenvalue corresponds to the 
homogeneous mode which is described by the coarse- 
grained operator A$ = X^ugs ^m/VpI- F° r 1^1 = O(Z), 



we obtain a macroscopic occupation Ns = (AgAg) 3> 1 
after a finite time t. The precise scaling depends on the 
size |<S|. For |<S| -C c 2 t 2 , the condensate fraction Ns/\S\ 
is basically independent of |<S| and grows with exp{ 7 t}. 
For |<S| 3> c 2 t 2 , on the other hand, the condensate frac- 
tion decays ~ This suggests that several frag- 
mented condensate patches of size 0(ct) form within the 
region S which are not yet fully coherent. 

These findings motivate the study of the phase corre- 
lations. To this end, we introduce a coarse-grained phase 



operator tps via A$ = exp{iips}V Ns with Ns = A^As- 
For a macroscopic occupation (Ns) ~> 1, the number 
fluctuations are negligible so that Ns ~ N$. Conse- 
quently, the phase between two non-overlapping regions 
S and S' of sizes 1 <C |<S|, |<S'| c 2 t 2 correlates according 



4 



to (for large 7*) 

f 7 Ar 2 1 

(exp{i(tp s ~ <fs')}) « exp j — ^ , (15) 

where Ar is the distance between the two regions. The 
distance over which the phase correlation spreads obeys 
a diffusion-like law Ar 2 ~ c 2 t/^f. For smaller distances, 
the relative phases become locked. 

As the final example, the four-point correlation 

(%ala R a\) = Tr^ KA {p^uK\alala K a\} = 
= (fit a K ) (at a x ) + (a],a A ) {a) v a K ) + 0(l/Z 3 ) , (16) 

is completely determined by (|T4l in leading order 1/Z 2 . 
Despite the similarities, this result is not based on the 
usual Wick expansion (we have a strongly interacting 
theory) but on the hierarchy © and the properties of 
the initial Mott state. For example, for (a K ) 7^ 0, we 
would get an additional contribution of the same order 
from p c ^ x = 0(1/Z 2 ). The above equation allows us to 
calculate the long-range current correlation 

For the sake of completeness, let us briefly discuss a 
quench deep into the superfluid regime J > U, where a 
finite wavenumber fc* yields the fastest growth. In this 
case, the correlations behave in a way which is drastically 
different from ((14)) . The temporal growth and the spatial 
dependence factorize, see also (3) 

(4(t)a M (i)) «M(t)exp{7**}.7o(fc.|TV-r,,|), (18) 

with the Bessel function Jq and the normalization factor 
Af*(t) = 0(1/ Z). As a result, the correlations grow but 
do not spread as before in (fl~4)) . i.e., with some velocity. 

In summary, we derived a hierarchy of correlations 
@ in order to describe the non-equilibrium dynamics 
of a lattice Bose gas ([5]) based on the expansion in 
inverse powers of the large coordination number 1/Z. 
The lowest order coincides with the Gutzwiller approach, 
cf. Eq. ©, and the higher orders describe the correla- 
tions, cf. Eq. (J7J. This method is applied to calculate 
the creation and amplification of quantum correlations 
in a quenched Mott-superfluid phase transition. We find 
that the off-diagonal long-range order (fT4]) spreads with a 
constant velocity and obeys universal scaling laws. The 
correlator (TT5l) of the phase associated with local con- 
densate patches expands with a diffusion-like law (phase 
locking). As an example for higher-order correlations, we 
calculated the four-point function (fl"6f which yields the 
current correlator (fl~7|). 

This work was supported by the SFB/TR 12 of the 
German Research Foundation (DFG). Helpful discus- 
sions with I. Cirac, F. Queisser, K. Krutitsky, and M. Pa- 
ter are gratefully acknowledged. 



[1] T. W.B. Kibble, J. Phys A 9, 1387 (1976); W.H. Zurek, 

Nature 317, 505 (1985). 
[2] A. Osterloh, L. Amico, G. Falci, R. Fazio, Nature 416, 

608 (2002); T. S. Cubitt and J. I. Cirac, Phys. Rev. Lett. 

100, 180406 (2008). 
[3] L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalattore, 

and D.M. Stamper-Kurn, Nature 443, 312 (2006); C.N. 

Weiler, T.W. Neely, D.R. Scherer, A. S. Bradley, M.J. 

Davis, and B. P. Anderson, Nature 455, 948 (2008). 
[4] M. Greiner et al, Nature 415, 39 (2002); F. Gerbier et 

al, Phys. Rev. A 72, 053606 (2005); W. S. Bakr et al, 

larXiv: 1006 .07541 
[5] E. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 

(2002); W.H. Zurek, U. Dorner, and P. Zoller, ibid. 

95, 105701 (2005); B. Damski, ibid. 95, 035701 (2005); 

J. Dziarmaga, ibid. 95, 245701 (2005); I. Klich, C. Lan- 

nert, and G. Refael, ibid. 99, 205303 (2007); B. Damski 

and W.H. Zurek, ibid. 99, 130402 (2007); A. Lamacraft, 

ibid. 98, 160404 (2007); D. Sen, K. Sengupta, and 

S. Mondal, ibid. 101, 016806 (2008); M. Moeckel and 

S. Kehrein, ibid. 100, 175702 (2008); F. M. Cucchietti, 

B. Damski, J. Dziarmaga, and W. H. Zurek, Phys. Rev. 

A 75, 023603 (2007); H. Saito, Y. Kawaguchi, and M. 

Ueda, ibid. 75, 013621 (2007); ibid. 76, 043613 (2007); 

A. Polkovnikov, S. Sachdev, and S.M. Girvin, ibid. 66, 

053607 (2002). 

[6] R. Schiitzhold, J. Low Temp. Phys. 153, 228 (2008); 
A. Polkovnikov and V. Gritsev, Nature Physics 4, 477 
(2008); A. Polkovnikov, Phys. Rev. B 72, 161201(R) 
(2005); J. Dziarmaga. larXiv : 0912 . 4034fr 2. 

[7] S. Sachdev, Quantum Phase Transitions (Cambridge 
University Press, Cambridge, England, 2001). 

[8] M. P. A. Fisher, P. B. Weichman, G. Grinstein and 
D. S. Fisher, Phys. Rev. B 40, 546 (1989). 

[9] R. Schiitzhold, M. Uhlmann, Y. Xu and U. R. Fis- 
cher, Phys. Rev. Lett. 97, 200601 (2006); U. R. Fischer, 
R. Schiitzhold, M. Uhlmann, Phys. Rev. A 77, 043615 
(2008). 

[10] R. Schiitzhold, M. Uhlmann, and U. R. Fischer, Phys. 

Rev. A 78, 033604 (2008). 
[11] See, e.g., M. Christandl, R. Koenig, G. Mitchison, 

R. Renner, Comm. Math. Phys., 273, 473 (2007); and 

references therein. 
[12] R. Kubo, J. Phys. Soc. Japan 17, 1100 (1962). 
[13] R. Balescu, Equilibrium and Nonequilibrium Statistical 

Mechanics (Wiley, New York, 1975). 
[14] K. Krutitsky, P. Navez. larXiv: 1004.2T21 , D. van Oosten, 

P. van der Straten, and H. T. C. Stoof, Phys. Rev. A 

63, 053601 (2001); D. van Oosten, D. B. M. Dicker- 

scheid, B. Farid, P. van der Straten, and H. T. C. Stoof, 

Phys. Rev. A 71, 021601(R) (2005). 
[15] D.S. Rokhsar and B.C. Kotliar, Phys. Rev. B 44, 10328 

(1991); K. Sheshadri, H. R. Krishnamurthy, R. Pandit, 

and T. V. Ramakrishnan, Europhys. Lett. 22, 257 (1993). 
[16] M. Uhlmann, R. Schiitzhold, and U. R. Fischer, Phys. 

Rev. Lett. 99, 120407 (2007); Phys. Rev. D 81, 025017 

(2010). 

[17] E. H. Lieb and D.W. Robinson, Commun. Math. Phys. 
28, 251 (1972). 

[18] L. P. Pitaevskii and S. Stringari, Bose-Einstein Conden- 
sation, Oxford University Press, Oxford, 2003. 



